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We consider the ground state of a harmonically trapped Bose-Einstein condensate within the 
Gross-Pitaevskii theory including the effective-range corrections for a two-body zero-range poten- 
tial. The resulting non-linear Schrodinger equation is solved analytically in the Thomas-Fermi 
approximation neglecting the kinetic energy term. We present results for the chemical potential and 
the condensate profiles, discuss boundary conditions, and compare to the usual Thomas-Fermi ap- 
proach. We discuss several ways to increase the influence of effective-range corrections in experiment 
with magnetically tunable interactions. The level of tuning required could be inside experimental 
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I. INTRODUCTION 

The Gross-Pitaevskii (GP) equation 0, @, Q has been 
extremely successful in describing a wide range of mean- 
field features for experiments with Bose-Einstein conden- 
sates (BECs). Inparticular, the Thomas-Fermi (TF) 
approximation [3, 0, [f| , where the kinetic energy is ne- 
glected, has been very rewarding 0] • This approximation 
holds for repulsive condensates with positive scattering 
length a and large particle numbers. In the regime of 
validity of the TF approximation, the total energy is dis- 
tributed between interaction energy and potential energy 
from the confining trap, while the kinetic energy becomes 
negligible. 

Because of the non-linear nature of the GP equation, 
it is only solved analytically in a few cases, e^., vortices 
and solitons in homogeneous condensates 2, 31. The TF 
solution is also analytical, although it only holds in the 
bulk of the condensate. At the surface the approxima- 
tion breaks down and is usually patched by including the 
kinetic energy at the surface [5j, [6| • 

The interactions of the ordinary GP equation are based 
on the lowest order zero-range potential, which is gov- 
erned by the scattering length alone. Although this ap- 
proximation is usually very good, the higher-order correc- 
tions to the scattering dynamics [1, [l(| can be crucial 
in certain cases, e.g., for Rydberg molecules embedded in 
BECs [lfj and for narrow Feshbach resonances In- 
clusion of higher-order terms is well known and applied in 
Skyrme-Hartree-Fock calculations in nuclear physics [l^] ■ 
Here, they often play a crucial role in order to get bulk 
nuclear properties right [HI, [3. However, the effects of 
similar higher-order terms in the GP equation have been 
less investigated. 

In this paper, we solve the modified GP equation with 
higher-order interactions analytically in the TF approxi- 
mation. The paper is organized as follows. In Sec. [Til wc 
introduce the modified GP equation and its parameters 
and show how it is derived from an appropriate energy 
density functional with careful treatment of boundary 
terms. We present the analytical solution in the TF ap- 



proximation in Sec. Hill and discuss the condensate size 
and chemical potential as function of the interaction pa- 
rameters in Sec. IIVI The density profiles and energies 
are discussed in Sec. [V] and in Sec. IVI1 we address the 
consistency of the TF approximation by considering the 
kinetic energy of the solutions. We compare to some rel- 
evant atomic systems in Sec. IVUl and finally present our 
conclusions in Sec. IVIIII 



II. MODIFIED GP EQUATION 

We assume that the condensate can be described by 
the GP equation. Since we are interested in the ultra- 
cold regime, where the temperature is much smaller than 
the critical temperature for condensation, we adopt the 
T = formalism. In order to include higher-order effects 
in the two-body scattering dynamics, we use the modi- 
fied GP equation derived in [lfj , which in the stationary 
form reads 
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where m is the atomic mass, V is the external trap, 
Uq = 4irh 2 a/m, and gi = a 2 /3 — ar e /2, with a and r e 
being, respectively, the s-wave scattering length and ef- 
fective range [Toj | . We assume an isotropic trap, V(r) = 
raw 2 r 2 /2, and introduce the trap length b = y^h/moj. 
The single-particle density, p(r) — \^(r)\ 2 , is normalized 
to the particle number, N = J drp(r), and \x is the chem- 
ical potential. 

As the boundary conditions are important for the TF 
approximation applied below we now discuss the proce- 
dure for obtaining the modified GP equation from the 
corresponding energy functional which is 



= J dr(e K + e v + £/ + e Z2 ), 
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The corresponding integrated energy contributions are 
denoted Ek, Ey, Ej, and Ej2, respectively. To obtain 
Eq. U), we vary Eq. ([2]) with respect to \I>* for fixed 
To first order in S^f*. we have 
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Here, S is the outward-pointing surface normal. In the 
usual analysis, one assumes that \& and vanishes at 
infinity, drops the boundary terms, and Eq. |T]) is ob- 
tained by varying E — pN . However, the existence of 
these surface terms is essential for the inclusion of higher- 
order interactions as discussed below. 

In the rest of this paper we use trap units, hw = b = 1, 
i.e., energies (E, V, p, etc.) are measured in units of Yuo 
and lengths (a, r e , r, etc.) in units of b. Note that 172 has 
dimension of length squared. 



III. THOMAS-FERMI APPROXIMATION 

Let us briefly review the standard Thomas-Fermi ap- 
proximation [J 0, H, El- Neglecting the kinetic-energy 
term, as compared to the trap and interaction energies, 
the GP equation has the solution 
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with chemical potential ptf- This solution is used out 
to the surface, Rtf, while outside ptf = 0. The nor- 
malization and surface condition ptf(Rtf) = give 



Ptf - -^Rtf- 

The total energy becomes 

Etf 
N ~~ 



Rtf = (15iVa) 1/5 
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The trap and interaction energies are Ey — 3E/5 and 
Ei = 2E/5, respectively. Since Rtf > in Eq. 0, 
these results only hold for a > 0. The TF approximation 
is good for Na ^> 1, except at the surface region where 
the kinetic-energy density diverges. Here, the solution 
can be corrected as in @, H, H, essentially giving a 
small exponential tail. 



We now consider the TF approximation with the 
higher-order interaction term, ej 2 . Ignoring the bound- 
ary terms in Eq. ^ , the modified GP equation can then 
be written in terms of the density p(r) = \^(r)\ 2 as 



fi = -r 2 + Ana (p + g 2 W 2 p) 
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With scaled coordinate x — r / \fg2~ (assuming g 2 > 
for the moment) and density f(r) = Anaxp(r)/ g 2 , this 
becomes 
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The inhomogeneous and homogeneous solutions with 
boundary condition /(0) = are 
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fi(x) = ( -x + 3)x, fh(x) = — sinx, (11) 
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where A is a constant (with dimensions of length squared) 
to be determined later. The full solution is 
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For a given A, the chemical potential p and the conden- 
sate radius R are determined by the normalization and 
the surface condition, 



iwx 2 p(x)dr — N and p{xq) = 0, 



(13) 



where £0 = Rl \f§2- The solution p should be positive 
for x < xo which must be explicitly checked. Outside xq, 
we use p = 0. 

We now consider the boundary terms in Eq. 
Above, we assumed that p(xo) = at some finite radius 
xq which we identify as the condensate size. However, 
only the first two boundary terms in Eq. ((5|) vanish on 
account of this condition. For the last term in Eq. © to 
vanish we need V x ^(xq) — 0, which implies that 



dp 

dx 



{x ) = 0. 
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Notice that this latter derivative is in fact non-zero in the 
g 2 = case, which is the root of the divergence of the 
kinetic energy at the condensate surface as we discuss 
later. Equation (Tj"4"]) gives a closed expression for the 
remaining free parameter A, 
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Xq COS Xq — Sill Xq 
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This additional requirement on the derivative at the 
edge of the condensate implies that higher-order terms 
require a smoothing at the surface of the cloud. In ad- 
dition, the discussion of which kinetic operator structure 
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to use (|V\1/| 2 or \1/*V 2 \1/ @) is obsolete in our treatment 
since the boundary term d4 f *V4 r vanishes. In this sense 
the inclusion of a higher-order term neatly removes some 
of the difficulties of the traditional TF treatment. 

The solutions with a finite boundary R of the modi- 
fied GP equation only minimize the energy functional if 
Eq. (TH| holds. We note that extremal states of the en- 
ergy functional always satisfy the virial theorem. Thus, 
enforcing the virial theorem on the GP solutions is equiv- 
alent to Eq. (TH| . We show in the Appendix that the 
virial theorem approach also leads to Eq. (fT5|) . 
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IV. SIZE AND CHEMICAL POTENTIAL 
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We now determine the condensate size R and chemical 
potential ji. The normalization condition is 



Na _ 3 / fi x% 
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while the surface condition reads 



xl/2 + 3 + 
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Combining eq. p5 |) - (fT7|) gives 
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which determines R for given Na and g 2 , and upon back- 
substitution also \i. 

The 52 < case can be worked out analogously by 
replacing trigonometric functions with hyperbolics and 
keeping track of signs. The two cases can in fact be 
combined into one equation 
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This equation determines Xq = R 2 / g 2 implicitly as 
function of Na/\g 2 | 5 ^ 2 . The result is shown in Fig. [T] We 
notice that in principle, R becomes a multi-valued func- 
tion. However, all the higher solutions for g 2 > [dotted 
in Fig. [T] are spurious, since the density becomes nega- 
tive on one or more intervals inside R. The non-spurious 
solutions [solid line in Fig. [T] define R as a single- valued 
function of a and g 2 , which was not guaranteed a priori. 
The four quadrants in Fig. [1] correspond to the different 
sign combinations of a and g 2 . The sign of the extra in- 
teraction energy, Ei 2 , is determined by ag 2 V 2 p. For a 
typical concave density, the Laplacian term will be nega- 
tive. We therefore see that for ag 2 > 0, the higher-order 
interaction is attractive, whereas for ag 2 < 0, it is repul- 
sive. The TF solution only exists for ag 2 < 0. We discuss 
both cases separately below. 



FIG. 1: (Color online) Condensate size (R) as function of Na 
and i?2 as found in the modified TF approximation, Eq. (|19|l . 
The solutions (a) and (b) correspond to the sign combinations 
(a > 0, cj2 < 0) and (a < 0, t/2 > 0), respectively. No 
solutions exist for ag2 > 0. The spurious solutions (dotted) 
have negative densities for one or more intervals inside R. 
The branch (a) approaches the normal TF result Eq. when 
Na — > +00 or (72 — » —0. Note that the convergence is only 
relative [see eq. (|19|> ] and the TF limit is better represented 
in the logarithmic inset. Points indicate the data from Tab. [I] 
All values are in trap units. 



A. The attractive regime: a<?2 > 

For a < 0, g 2 < [third quadrant in Fig. [I] there are 
no solutions, which is expected since the normal TF ap- 
proximation has no solutions for a < as the interaction 
energy Ej is negative and the kinetic energy that could 
prevent collapse is neglected. 

The g 2 > 0, a > case in the first quadrant has 
only spurious solutions. Here the g 2 term is attractive 
for the typical concave density and a collapse towards 
a high-density state is possible in complete analogy to 
the usual discussion of attractively interacting conden- 
sates within the standard GP theory. Whereas there can 
be metastable states at large values of Na/g 2 ^ 2 , these 
are stabilized by kinetic energy and thus are not present 
in our TF approach. Thus, even when the total kinetic 
energy is small, it is still needed to prevent the attrac- 
tive higher-order term from amplifying local-density vari- 
ations. 

This important point can also be established by con- 
sidering the stability of the homogeneous condensate 
through linearization of the GP equation. By repeating 
the analysis of [2| with the higher-order term, we find 
that for g 2 > and a > 0, the kinetic-energy term is 
crucial for the stability of the excitation modes. In fact, 
exponentially growing modes will always be present if the 
kinetic energy is neglected. This will be discussed else- 
where in relation to the numerical solution of the full GP 
equation (l5| . 
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Na/\g 2 \ 5 / 2 

FIG. 2: (Color online) Chemical potential g, as function of 
Na and <?2 as found in the modified TF approximation, using 
the solutions (a) and (b) from Fig. [T] For branch (a) and the 
upper part of branch (b) (see inset), we have g > 0. The 
lower part of (b) has g < 0. Points indicate data from Tab. [J 
All values are in trap units. 



B. The repulsive regime: agi < 

For <72 < 0, a > a single solution (a) exists. This 
was expected since £72 > gives extra stability. The so- 
lution approaches the normal TF result in Eq. when 
iVa / 1 <72 1 5 / 2 — > +oo, as can also be seen from Eq. (fTO)) . Of 
course in this limit E12 <C Ej. However, the —1 term in 
Eq. (JTHJ) implies that the convergence to the normal TF 
solution is only on a relative scale and is better repre- 
sented on a logarithmic scale as in the inset in Fig. [1] 

For gi > 0, a < there is a single solution (b) 
which connects smoothly to the (a) solution. In the limit 
Na/\g2\ 5 / 2 — > — oo, which is determined by xo cotx = 1, 
we find R 2 /g2 = 20.1907. This solution is possible when 
the g% term provides just enough repulsion to cancel the 
usual a < collapse behavior. 



C. Chemical potential 

In Fig. [5] we show the chemical potential for the 
smoothly connecting solutions (a) and (b). Again we 
see that (a) approaches the normal TF limit for large 
Na/\g2\ 5 ^ 2 ■ Here, it is interesting to note how /i turns 
around near the origin [amplified in the inset in Fig. [5] 
and maintains a positive value. This occurs in the region 
where the lowest-order interaction gives a large negative- 
energy contribution which the 172 term is still able to bal- 
ance yielding a well-defined TF solution. This behavior 
is analogous to the balancing of attraction by the kinetic 
term in the usual a < 0,g 2 — case 0, Q. As a be- 
comes increasingly negative, so too does fj, and collapse 
is inevitable (and likewise when g2 — > + ). 
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FIG. 3: (Color online) Densities for branch (a) in Fig. [T] (172 < 
and Na = 10 4 ). The <?2 = —0.1 curve is on top of the normal 
TF result. Inset shows the smooth behavior at the surface for 
gi < 0. All values are in trap units. 

V. DENSITIES AND ENERGIES 

With R and g, determined, we can find the density pro- 
file, energy densities and integrated energy contributions. 
With Eq. (fl~2j) . the energy densities are given by 

ev = ej = 27rap 2 , (20) 

ei2 = -kp(3 + £ S -T)- (21) 

Using Eq. ((9]), the total energy density (without ck) be- 
comes 

e = ev + ei + e n = ~p(x)(y(x) + £). (22) 
2 32 

In Fig. [3l we show the density profile of the (a) solu- 
tions for Na — 10 4 and selected 32 < 0. We clearly see 
that the higher-order term tends to expand the conden- 
sate through its repulsion. Importantly, at the boundary, 
there is a smoothing caused by the condition in Eq. (fl4|) 
[see inset in Fig. [3]. We will discuss how this affects 
the estimated kinetic energy in the next section. As I32I 
grows, we see the condensate flatten and in the limit of 
very large I32I, it becomes a constant density. 

Figure 2] displays the density profile for the (b) solu- 
tions with a < for selected 32 > 0. Here, we see the 
profile collapse toward the expected delta-function with 
decreasing g 2 . It is interesting to follow the (a) solution 
through the origin in Fig. [T] and onto solution branch (b), 
passing from g 2 — —00 to 172 = 00. On the (a) branch, 
the solution flattens as g 2 decreases and eventually be- 
comes effectively constant in space. This is also true for 
the (b) branch at 172 = 00, and as 32 is decreased, the 
solution proceeds to shrink as the 32 term becomes un- 
able to provide the repulsion needed to prevent the a < 
collapse induced by the lowest-order term. 

From the figures, we see that large l^l induces large 
changes in cloud size. As the condensate can be imaged 
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FIG. 4: (Color online) Same as Fig. [T] but for solutions (b), 
i.e., opposite signs gi < and Na = — 10 4 . 



92 > 



92 < 







TF-limit a- 




K' 


N\a\ = 10 4 


E v /N 

E r /N 

En/N - 

ltfxE K /\E\ ■■ 
i i i i 



40 
30 
20 
10 

-10 
-20 
-30 

-300 -200 -100 100 200 300 400 500 

Na/\g 2 \^ 

FIG. 5: (Color online) Different total-energy contributions. 
N\a\ = 10 4 . All values are in trap units. 



with very good resolution Q , this should be measurable 
if the regime of large I32I can be accessed. 

We now discuss the energy contributions which are in- 
teresting since the release energies are in fact measur- 
able quantities 0]. Since we neglect the kinetic term 
in the TF approximation, the release energy is simply 
Er = Ei + E12 = E - E v . In Tab. HI we give the in- 
tegrated energy contributions for some relevant values 
of 92 calculated for N\a\ — 10 4 , whereas Fig. [5] gives 
the energies as function of Na/\g2\ 5 ^ 2 ■ We note that for 
smaller values of 7V|a|, the same overall behavior is found, 
however, the kinetic term is more important and the TF 
approximation becomes worse. 

We observe that E/N grows towards the \g%\ = 00 
point. This is due to the trap energy increasing as the 
density flattens [Ey diverges around the origin in Fig. [5] . 
Furthermore, as 92 — > + the energy diverges toward —00 
as the collapse sets in [Ej diverges on the 92 > side 
in Fig. [5] . The boundary where the energy vanishes is 



around 32 = 5.14 for N\a\ = 10 4 , but this depends on 
the choice of N\a\. With respect to the release energy, 
we find that somewhere in the region 10 < 92 < 50, 
En becomes negative. This is a result of the unavoidable 
collapse, and also indicates that kinetic energy cannot be 
ignored at this point. Notice, however, that the release 
energy changes considerably and could provide a way to 
measure the influence of the 92 term. 



VI. CONSISTENCY OF THE THOMAS-FERMI 
APPROXIMATION 

We now address the validity of the TF approximation 
with the 32 term included. In order to do so, we must 
consider the contribution of the kinetic energy. The ki- 
netic energy density can be written as 



92 
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A x sin x — cos x 
92 x 



(23) 



Strictly speaking, this is not the true kinetic energy, since 
the kinetic terms were neglected from the start. However 
Eqs. and can be used to test whether the TF 
approximation holds locally, i.e., <C e should hold for 
the solution p to be consistent. In Tab. Q] we calculate 
the integrated contribution of the kinetic energy relative 
to the total TF energy and we find that the contribution 
is small everywhere except the point where E = on the 
92 > side of Fig. [5] Here, the kinetic energy is of course 
the most important term and the TF approximation is 
poor. 

In the standard TF, the kinetic energy causes trou- 
ble at the boundary of the cloud. Here, V'F cx Vp/^/p 
and since the density vanishes and the derivative is finite 
[see Eq. [6], this diverges at Rtf- When including the 
higher-order term we need to use the additional bound- 
ary condition = at R, so the kinetic energy will 
be strictly zero at R. However, as one approaches the 
boundary, the kinetic- energy density grows rapidly be- 
fore it descends towards zero within a very small interval 
at R. The total energy density in Eq. ([22]) goes to zero 
at this point and we find that ex /e is very large near the 
boundary as in the usual §2 = case. 

We conclude that the inclusion of the higher-order 
term does not alleviate the difficulties with kinetic en- 
ergy at the boundary. The techniques for addressing this 
problem described in 0, Q should therefore be general- 
ized to include the higher-order interaction term in order 
to improve the description at the boundary of the cloud. 



VII. COMPARISON TO ATOMIC SYSTEMS 

The considerations above show that deviations from 
the usual TF approximation can be strong when 172 is 
large. In the following, we reintroduce explicit units for 
comparison with real systems. We have to consider 92 /b 2 . 
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"The kinetic energy estimated by surface corrections as in [2j. 
''Values are indicated by points in Figs, [l] and [2] 
c The total energy \E\ is zero near 32 = 5.14, hence the TF ap- 
proximation is invalid here. 



TABLE I: Chemical potential n and condensate size R for different and fixed 7V|a| = 10 4 . The integrated energies are trap 
(Ey), interaction (Ei, E12), total [E = Ey + Ei + E12), and release energy (Er — E — Ev)- The TF limit is approached for 
Q2 — * —0. The ratio of kinetic energy Ek to total energy E indicates where the TF approximation is valid. The corresponding 
density distributions are shown in Figs. [3] and [4] All values are in trap units. 



Of course, the b 2 factor means that this quantity is gen- 
erally very small since g-2. is of order Oq and b is of order 
10 4 a . 

We first consider some typical background values for 
bosonic alkali atoms away from resonance. We estimate 
the effective range to be the of order of the potential 
range and assuming a van der Waals interaction, we have 
r e ~ 50 — 200arj. For typical one-component gases we 
have -450a < a < 2500a 0. Since g 2 = a 2 /3-ar e /2, 
we see that the a 2 term will dominate and in all cases 
< 92 ^ 10 6 a . In trap units, this becomes g 2 /b 2 < 
5 • 10~ 3 (l//m/6) 2 . In typical traps of b ~ 1 — 10/xm, the 
higher-order term is therefore very small. These values 
also predominantly lie in the first quadrant of Fig. [T] and 
thus no TF solution exists. 

Let us first consider Feshbach resonances in order to 
increase the influence of the g 2 term. We use a multi- 
channel Feshbach model [l^l, which describes the full 
T matrix as a function of resonance position Bq , width 
A£>, magnetic-moment difference between the channels 
A/i, and the background scattering length a& g . Per- 
forming an effective-range expansion (TlT |. we have a = 
a bg [l-AB/(B-B )] and r e = r e0 /[l - (B - B )/AB} 2 , 
where r e o = —2h 2 /mai,gAiiAB < 0. Combining these 
relations, we find r e = r e o(l — a bg/o) 2 and 



a ar e0 , db g , 2 

52(a) = — — (1 a -) . 

6 2 a 



(24) 



Hence g 2 diverges when a — > (referred to as zero- 
crossing) or a — > 00 (on resonance). Near zero cross- 
ing, the effective-range expansion is, however, severely 
divergent and its validity is questionable. Even so, the 
effective-range corrections near zero-crossing obtained 
are in fact identical to those obtained from use of the 
full T-matrix 18]. One finds lim a ^o a <?2 = keo|a| ff /2, 
where r e o < 0. 



As a concrete example, we consider the alkali isotope 
39 K where several Feshbach resonances of vastly different 
widths were found recently [l!|. First, we focus on zero- 
crossing and consider the very narrow resonance at Bq = 
28.85G with AB = -0.47G, A/i = 1.5/x B , and a bg = 
— 33ao- We obtain r e0 = — 5687ao and 052 — ► 93.8 • lO 3 ^ 
for a — > 0. It is important to notice that 0172 > around 
a = 0. This means that we are looking for solutions 
in the first and third quadrants of Fig. Q] and again we 
have to conclude that no TF solutions can be found when 
higher-order terms are taken into account. 

Another case of interest is around resonance where 
\a\ = 00. Here, we have r e ~ r e Q and g 2 cx a 2 > 
on both sides of the resonance. Thus, the a > side 
will be in the first and the a < in the second quadrant 
of Fig. [T] This makes it difficult to imagine sweeping 
the resonance from either side to probe the solutions on 
branch (b) in Fig. [TJ One could imagine starting on the 
a > side with small g 2 > 0. The full GP equation will 
have perfectly sensible solution here, however, when one 
approaches the resonance the g 2 term will diverge and in- 
duce collapse already on the a > side. If we approach 
from the a < side, then we face the problem that the 
critical number of particles decreases dramatically before 
g 2 grows sufficiently and one therefore needs a very small 
condensate since Na/b ~ 0.5 [ll|. At this point, the TF 
approximation is no longer valid. 

The Feshbach resonance used to increase 172 must be 
very narrow in order for r e o to be large. However, most 
experimentally known resonances are not narrow. For 
broad or intermediate resonances, we have to consider 
the long-range van der Waals interaction when calculat- 
ing the effective-range corrections. Analytic formulas for 
this case have been worked out in [2(| , and we note that 
the effective range diverges as a~ 2 near zero crossing ex- 
actly as in the Feshbach model above. For very narrow 
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resonances, we still have <C r e o, where j3§ is the char- 
acteristic length of the van der Waals interaction. The 
model above should thus give the dominant contribution. 

Using the van der Waals formulas we can estimate agi 
at zero crossing. We find 

where x e — (T[l/4]) 2 /27r, with T the gamma function. 
We have explicitly introduced the oscillator length which 
is the relevant length scale of comparison. Importantly, 
we find that agi < for a > and we are thus in the 
fourth quadrant where a TF solution exists. For a < 0, 
we pass to the second quadrant as g 2 becomes positive 
and a single collapsed solution can be found. 

We now estimate the parameters obtained from the van 
der Waals formulas. With b = 1/im and (3q ~ 123. 3ao Q, 
we have ag 2 /b 3 ~ — 1CU 8 . We thus have Na/\g 2 \ 5 ^ 2 oc 
10 8 (iVa)a 5 / 2 . For values of a that are not extremely 
small, the solution is therefore typically located far to 
the right in Fig. [1] where it will look similar to the g 2 = 
case. We can estimate how close to zero one would have 
to tunc a in order to see deviations using the a — > 
limit of the van der Waals effective range. Let us aim for 
92 /b 2 = —10 which should be observable in the conden- 
sate profile according to Fig. [3] With b = 1/im, we need 
a ~ 10~ 6 /?6 ~ 1.7 x 10~ 4 ao. Using broad resonances, 
one can tune to zero at the level of 10 _2 ao in 39 K [2l| . 
Observing the effect of the g2 term therefore seems out of 
reach at the moment, but might be possible in the near 
future. Of course, we still have to maintain a large value 
of Na for kinetic energy to be small, and thus a larger 
condensate is needed close to zero crossing. 

From the examples above, we see problems in accessing 
the TF solutions presented above in current experiments 
with ultracold alkali gases. In particular, we notice that 
realistic systems which have been used for creation of 
BECs in alkali-metal gases for the last decades have pa- 
rameters that predominantly lie in the first quadrant of 
Fig. [TJ As we have discussed, there are no well-defined 
TF solutions in that region. Therefore, we see that the 
kinetic energy plays a decisive role and we are forced to 
consider it in principle, even if it is small for all practi- 
cal purposes. The physical reason is that for a > and 
g 2 > 0, the higher-order interaction is effectively attrac- 
tive and induces collapse which will have to be balanced 
by a barrier from the kinetic term, similar to the a < 0, 
92 = case [l|. Since we neglect the kinetic term in the 
TF approximation, we should not expect to find solutions 
in the ag2 > case. 

Only in the case of resonances dominated by the long- 
range van der Waals interaction do the parameters allow 
for TF solutions with non-zero g 2 . However, here the 
length scale of the trap makes the contribution very small 
and the TF solution becomes identical to the g 2 case. 
One could in principle tune a very close to zero-crossing 
and obtain a significant contribution but the level of tun- 
ing required is beyond current experimental reach. 



VIII. CONCLUSIONS 

We have considered the effect of higher-order inter- 
actions in Bose-Einstein condensates within the Gross- 
Pitaevskii theory. We derived the GP equation with 
effective-range corrections included and solved it analyt- 
ically in the Thomas-Fermi approximation. Higher-order 
interaction terms act as derivatives on the condensate 
wave function which means that the boundary conditions 
on the solutions of the GP equation must be carefully 
considered. We then discussed the solutions for various 
parameters and presented the chemical potential, density 
profiles, and the energy contributions. 

We find that no TF solutions are possible when the 
higher-order term is attractive. This conclusion holds 
both in the trapped system and in the homogeneous case 
(TH . An estimate of the relevant parameters for alkali 
atoms showed that away from resonances, they typically 
lie in the region where the effective-range correction is 
effectively attractive and likewise near very narrow Fesh- 
bach resonances. We conclude that in those cases, the ki- 
netic energy, even if very small, is crucial in order to sta- 
bilize collapse due to higher-order interaction terms. For 
broader resonances where the long-range van der Waals 
potential is dominant, we find that modified TF solu- 
tions exits. However, for typical traps, the parameters 
are very small and tuning of the scattering length near 
zero crossing at a level beyond current experimental reach 
is necessary. This might of course become possible as ex- 
perimental control improves in the future. 
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APPENDIX A: DETERMINATION OF A FROM 
THE VIRIAL EQUATION 

Even though Eq. (TT2"]) is a solution to the modified 
GP equation Eq. © for all A, it does not necessarily 
minimize the energy functional as discussed in the main 
text. This can also be seen from the virial equation (with 
neglected kinetic energy), 

- 2E V + 3£/ + 5E I2 = 0, (Al) 

which holds for all extremal points of the energy func- 
tional. Equation (|A1[) is derived from the energy func- 
tional using scaling arguments as in Q. 

As an example, consider the A = solution in Eq. (|12p . 
This solution has a chemical potential shifted by 3g 2 
compared to the 52 = TF result. But the density 
is unchanged and so is Ey and E[. Hence, the usual 
virial equation — 2Ey + ?>Ej = for g 2 — also holds 
for g 2 7^ 0. Since Ej 2 = — 3g 2 /2 7^ 0, the virial equa- 
tion Eq. (|A1|) is not fulfilled, and hence the A = solu- 
tion is not extremal. Below, we use the virial equation 
to calculate the value of A that minimizes the energy 
functional and the corresponding R and fx. We will also 
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prove that this condition is in fact equivalent to assume 
p(xo) = W x p(xo) = at the boundary. 

The general results for A, R, and p can be derived 
using the normalization and surface conditions Eq. (|13[) . 
and the virial equation Eq. (|A1[) . For convenience, we 
introduce the variables p — p/(3g2) + 1, A = A/(3g2), 

5 /2 

and c = Na/g 2 ■ The different energy contributions are 

£V = 3s / dxx 4 (/2 + A ), 

Jo 6 x 

Ej = 9s I" dxx 2 (p - — + A S -^) 2 , 
Jo 6 x 

f x ° ,9/ x 2 T sinir, , T sina;. 
Ei 2 = -9s / dxx 2 (fi - — + A 1 + A , 

OX X 

(A2) 



where s — g^ 2 /(2a). Direct integration of Eq. (IA2[) . 



insertion of p from Eq. (fl7|) , and some algebra gives the 
virial equation 



= -2E V + 3J5j + 5E I2 
= — ^— (Xq — 3^4(o;o cos xq — sin xq)Y 

Xq 



(A3) 



We immediately see that this is in fact equivalent to 
Eq. (JT5J). Therefore, the solution we have explicitly found 
above minimizes the energy functional with boundary 
conditions p(xo) = V x p(xo) = 0. More generally, when 
we solved the modified GP equation without consider- 
ing the boundary terms in Sec. we found a one- 
parameter family of solutions (parametrized by A). The 
virial theorem is merely a constraint on A for obtaining 
a minimum of E. 
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